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We examine the dynamics of a dissipative photonic Josephson junction formed by the weak cou- 
pling of an optical soliton in a nonlinear dielectric waveguide and a co-propagating surface plasmon 
along a parallel metal surface with a linear dielectric spacer. We employ a heuristic model with 
a coupling function that depends on the soliton amplitude, and consider two phenomenological 
dissipation mechanisms separately: angular velocity dissipation and population imbalance dissipa- 
tion. In the former dissipation mechanism, the system exhibits phase-slip phenomenon where the 
odd-7r phase modes decay into even-7r phase modes. The latter damping mechanism sculptures the 
phase-space significantly by introducing complex features, among which Hopf type bifurcations are 
notable. We show that some of the bifurcation points expand to stable limit cycles for certain 
regimes of the model parameters. 
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I. INTRODUCTION 



As the research on the interaction of light with metal- 
lic structures matures as a well-established technology 
called plasmonics pJE] , the coupling of surface plasmons 
to different light sources are being investigated, with the 
motivation that controlling surface plasmons offer the 
potential for developing different types of SP-integrated 
nanophotonic devices [4]-[6]. In particular, the coupling 
between SPs and confined light modes are widely inves- 
tigated. A recent proposal is based on the resonant in- 
teraction between the SPs on a metal surface and co- 
propagating soliton in a nonlinear dielectric medium [7J. 
Under a classical formulation, this system exhibits rich 
nonlinear dynamical features where the interaction de- 
pends on the soliton amplitude, as such it may be uti- 
lized to manipulate the SP propagation. As an addi- 
tional feature, it has been shown [5] that the system is 
akin to the bosonic Josephson-j unction of Bose-Einstein 
condensates [8HT2] so that similar and different nonlin- 
ear Josephson junction features may be realized in this 
optic-plasmonic system. 

While this was an exciting analogy, surface plasmons are 
subject to strong dissipative effects in the host metal and 
perfect Josephson junction is a too idealized model. In 
the present work, we aim to make a more faithful repre- 
sentation of the system by taking into account the dis- 
sipation. Remarkably, we find that the dissipation can 
bring some benefits if it could be introduced under con- 
trol. The junction can exhibit stable dynamics in spite of 
dissipation. Based on this motivation, we investigate here 
the dissipation effects of the coupled soliton-SP system 
based on the model introduced in Ref. [7J and then for- 
mulated as a nonlinear Josephson junction [8], We note 
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that dissipation effects in the Bosonic Josephson Junction 
(BJJ) were previously studied, [Hj which might be useful 
to refer for comparison. We employ a standard dynamical 
analysis of the system in the phase space representation, 
and investigate how different dissipation mechanisms and 
other model parameters affect the phase-space landscape. 
Bifurcations in the system are highlighted which may be 
induced in certain regimes of these parameters. 

II. MODEL AND THEORETICAL 
FORMULATION 

Construction of a classical theoretical model for a cou- 
pled soliton surface-plasmon system has been addressed 
recently in a number of publications. An elegant initia- 
tive was the heuristic model introduced by Bliokh and 
collaborators, [7J in which the interaction between the 
soliton and the surface-plasmon is formulated as a cou- 
pled non-linear/linear oscillator system, with the cou- 
pling parameter being dependent on the soliton ampli- 
tude. While being constrained by a number of simplify- 
ing assumptions, the model predicts stationary coupled 
modes of the soliton and the surface-plasmons under fea- 
sible parameter regimes. In a following work, station- 
ary and quasi-stationary solutions for the so-called soli- 
plasmons are investigated under an asymmetric coupling 
model [13 . In this paper, we will employ the heuristic 
coupled-oscillator model, albeit with an improved cou- 
pling parameter adapted from the recent formulations. 
Our main motivation is to incorporate the correct be- 
havior for the coupling in the strong soliton and strong 
surface-plasmon amplitudes, respectively [13]. We be- 
gin by recapitulating the heuristic theoretical model in 
Ref. [7J, where the interaction between optical solitons 
and surface plasmons (SP) is discussed in the context of 
an resonant optical system: a soliton propagating par- 
allel to a metal interface can excite surface-plasmons by 
its evanescent lateral tail, which, in turn, interacts in 
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resonance with the soliton along the propagation. The 
soliton and the SP are assumed to be co-propagating, 
hence only the spatial dynamics of the propagation is 
of concern. The nonlincarity of the system is assumed 
to be confined at a distance, d, from the metal interface 
such that the surface-plasmon propagation retains linear- 
ity, and the coupling between the soliton and SP-fields 
is weak. The formulation is based on a 2D coordinate 
system in which y is the propagation direction and x is 
the lateral direction. 

The total electric field of the system is introduced by 
the following variational ansatz in which the soliton and 
SP fields are written in a product form of their respective 
longitudinal {c PtS (y)) and transverse amplitudes. 



E(x,y) =c p {y)e K * x + 



Evanescent wavevectors are k t 
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ky/j/2\c s \. The k is the incoming wave vector whereas k p 
is the SP wave vector. 7 is the nonlinearity parameter of 
the dielectric strip. The longitudinal amplitudes c PiS (y) 
obey the coupled spatial propagation wave equations: 



c p + /3 2 c p = q(\c s \)c s , c s + fi 2 s c s = q(\c s \)c p . 



(2) 



where c PjS = d Q^' s and £ = ky is the dimensionless prop- 
agation coordinate. (3 p — k p /k and /3 S = 1 + 7|c s | 2 /4 are 
the propagation constants of SP and soliton respectively. 
q{\ c s |) is the coupling function which we will discuss fur- 
ther below. The coupled equations are linearized by in- 
troducing c p s = C p s e^ into Eq. [2] and by employing 
the slowly varying amplitude approximation to discard 
higher order derivatives. The final set of equations are 
then: 
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v p = (3 p — 1 <C 1, v a = P s — 1 -C 1 are the small de- 
viations of the dimensionless propagation constants of 
surface-plasmon and soliton, respectively. 

Evidently, the coupling function plays the main role 
in these equations. Qualitatively, it should describe the 
interaction through the overlap of the evanescent tails of 
the soliton and plasmon in the (linear) dielectric spacer 
between the metal interface and the nonlinear dielectric 
waveguide. Since the lateral profile of a soliton depends 
on its amplitude, the coupling function is expected to ex- 
hibit this dependence, too. An analytical formulation of 
the coupling function was provided in Ref. [13] with the 
remark that the equation system is not symmetric in the 
coupling function in general. Within the heuristic model, 



the following functional form captures the expected be- 
havior in the strong soliton and strong SP amplitude lim- 
its: 



g(|C.|) * \C s \e 



rkdy/Z\C t \ 



(5) 



Here, a is an overall constant which we introduce as a 
parameter to account for the other constants (e.g. per- 
mittivities, propagation constants). The key parameters 
are preserved in the functional form. In the strong soliton 
amplitudes (\C S \ ~ 1), the exponentially decaying term 
dominates the coupling. In the weak soliton amplitude, 
the coupling is proportional to \C S \ since the soliton lat- 
eral profile in Eq. [T] would be wider resulting in a larger 
overlap. 

The analogy between this model and the bosonic 
Josephson junction (BJJ) dynamics [TOl H2] is revealed 
when we further substitute C SjP = C S}P e l ^ s - p into Eq. [4] 
and introduce a new variable set consisting of the frac- 
tional population imbalance Z — (|C S | 2 — \C p \ 2 )/N and 
the relative phase difference between the soliton and the 
SP 6 = 6 S - 6„; 



Z = -q(Z) y/l - Z 2 sin< 



AZ + AE- 



q(Z)Z 

VT^z 2 



(6) 
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where A = is the nonlinearity parameter, and AE = 
A — v p parametrizes the asymmetry between the soliton 
and SP states (similar to the asymmetry between the 



wells of a double well system). N = (\C S 



constant for the isolated system. The coupling function 
q takes the form 



9{Z) 



(1 + Z) crkdy / 2 A(i+z) 



(8) 



We stress that the Z dependence makes an inherently dy- 
namic coupling as opposed to constant (or externally tun- 
able) coupling parameter present in BJJ-systems. With 
this formulation, the coupled soliton SP system can be 
described as a nonlinear Josephson junction. 

The range of the model parameters are chosen to be 
around the resonant coupling regime [7]. In the calcula- 
tions presented below, we fix the dimensionless separa- 
tion parameter kd = 6, and choose the other parameters 
between A ~ 0.003 - 0.03, AE ~-0.4-0.4, accordingly. 
Evidently, using larger kd values would confine the inter- 
action to smaller soliton amplitudes and eventually result 
in the decoupling of the soliton and the SP completely. 
Small kd values would not be applicable in the weak cou- 
pling approximation. 



III. ANGULAR- VELOCITY DEPENDENT 
DISSIPATION 

Owing to the analogy to the BJJ models [5], the first 
dissipation mechanism we consider is about the incoher- 
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FIG. 1: (Colour online) (a)The Z — phase-space of the dissipationlcss SP-soliton coupled system with kd — 6, A — 0.01, AE — —0.0025. The 
magnitude of the gradient plot is color-coded logarithmically. White curves show typical open and closed trajectories, (b)Phasc-spacc trajectory 
with parameters of (a) and with angular velocity dependent dissipation, r\ — 0.2. Phase slip occurs between the cf> — — 7r and — modes. 



ent exchange of photons between the soliton and the SP. 
This is introduced by the term r]<fi in the following equa- 
tions, 

Z = -q{Z) y/l - Z 2 sin <j) - r)4>, (9) 
d> = AZ + AE + 4i£H= cos <b, (10) 

Vi-z 2 

Under coordinate reversal £ — > — £, we have Z — > Z 
and 4> ~~ ^ —4>i hence, the dissipative term rjcf) appearing in 
the equation for Z represents an irreversable process cor- 
rectly. The phase-space of the dissipationless and dissipa- 
tive system are shown in Fig[TJa) and[ljb), respectively. 
The gradient distribution is given with color-coded mag- 
nitudes and typical trajectories are indicated by white 
curves. The phase space is characterized by the zero- 
(even) phase and 7r-(odd) phase modes, where the criti- 
cal points are located. There is always one stable point 
at (j> = 0. At <p — i") one (stable) to three (two stable, 
one saddle) critical points may emerge depending on the 
values of the model parameters A,kd, and AE [5]. In 
Figure [l| a), sample trajectories (white curves) show an- 
harmonic closed orbits around stable critical points and 
open trajectories. The color indicates that the gradients 
get steeper in the SP-dominant population imbalance re- 
gion (Z < 0). 

Since the critical points are determined by the Z = 
0, <fi = condition, the </>-dependent dissipation does not 
alter the location of the critical points in the phase space 
but the eigenvalues of the Jacobian acquire nonzero real 
parts such that the stable critical points at <f> = it be- 
comes source and the critical point at <p = becomes a 
sink. While the trajectories flow towards the sink points, 
the phase-slip phenomenon is observed from </> =odd- 
7r-modes to <j> =even-7r modes, as reported for the BJJ 
systems [9]. Depending on the value of the asymme- 
try parameter, AE, the phase-slip occurs with a nonzero 
change in the averaged population imbalance < Z >. 



With a larger (smaller) dissipation constant, the typ- 
ical phase-diagram depicted in Fig[ljb) gets compressed 
(extended) along the phase axis but the features de- 
scribed above remain the same. Thus, while the tran- 
sient behavior depends on rj, the final state of the system 
is independent: ip s = <j) p and a fixed value of Z . 

IV. POPULATION IMBALANCE DISSIPATION 

In the following, we introduce the set of equations that 
incorporate a dissipation proportional to the population 
imbalance, Z. In the context of BJJ systems, such a dis- 
sipation arises from the finite lifetime of excited states of 
two condensates in different hyperfine levels in a single 
harmonic trap connected by tunnelling transitions |14j . 
In the analogous mechanical system of momentum short- 
ened pendulum, this corresponds to a damping propor- 
tional to the angular momentum variable [15] . which we 
consider to be applicable to the soliton-SP system phe- 
nomcnologically. 

Z = -q(Z) \fl - Z 2 sin - (Z, (11) 
A,= AZ+ AK+ q<yZ ^ Z cosd) (12) 

Vi - z 2 

The first thing we note here is that this dissipation 
term can modify the critical points in the phase-space 
significantly: their existence, location and corresponding 
Jacobian eigenvalues. For determining the critical points, 
using Z — 0, 4> = and the identity sin 2 (0) + cos 2 {4>) — 1 
we obtain the following root equation in Z , 

Z 4 C 2 + (l-Z 2 )[(AE+AZ) 2 {l-Z 2 )-q{Z) 2 Z 2 ] =0 (13) 

Figure[5]is a plot this equation for fixed kd = 6, A = 0.01 
and for (a) AE = 0, and (b) AE = -0.0025, with several 
order of magnitudes of the dissipation constant Blue 
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FIG. 2: (Colour online) (a)Plot of Eq.[l3]for kd = 6, A = 0.01 and AE = at various values of the dissipation constant Roots of f{Z) arc the 
critical points in the phase-space, (b)Plot of Eg. |13| for kd — 6, A — 0.01 and AE — —0.0025 at various values of the dissipation constant £. Roots 
of f(Z) arc the critical points in the phase-space. 



curve is plotted for dissipationless case (( = 0) for refer- 
ence. For AE = 0, Z = is always a critical point as all 
curves in Fig. ga) are tangent to Z = at this point, 
irrespective of dissipation. This can also be inferred from 
Eq. 13 Critical points close to Z m 1 are most sensitive 



to the dissipation strength and they already disappear at 
weak dissipation (( 1(T 4 ). For AE = -0.0025, as the 
dissipation gets stronger, the critical points collide and 
annihilate in pairs. This is actually indicating the pres- 
ence of saddle-node type bifurcations in the system which 
we show further below. Finally, above ( = 0.0275 no crit- 
ical point is present in the phase space. The sensitivity 
of the critical points to the dissipation strength increases 
rapidly after £ > 0.01. Table |T] classifies the critical points 
at different dissipation regimes for AE = —0.0025. Com- 
plex eigenvalues of the Jacobian with positive (negative) 
real part define spiral source (sink) points and saddle 
have zero eigenvalues. 



TABLE I: Classification of fixed points at different dissipation 
regimes for Fig. |2jb) 



C I 


II 


III 


IV V VI 


0<C<8.5 x 10" 5 saddle 


sso 


ssi 


saddle ssi saddle 


8.5 x 10" 5 <C<0.016 saddle 


sso 


ssi 


saddle - 


0.016<C<0.027 saddle 


sso 






0.027<C 









sso: spiral source, ssiispiral sink 

According to the typical ranges of the dissipation pa- 
rameter we now investigate the behavior of the system 
in the phase-space. Figure ga) is plotted for kd = 6, 
A = 0.01, AE = 0, and C = 10~ 2 . These values corre- 
sponds to the cyan curve in Fig. gb) . Z = is a spiral 
sink at <f> = and a saddle point at <p = n. Another spi- 
ral sink is located at Z = —0.28, 4> — 0.857T. The saddle 
point at Z — —0.98, <j> = 0.5ir is not visible at this scale. 

In Figure gb), same parameters as in Fig. 3^sl) are 
used but AE — —0.0025. Remarkably, a stable limit 



cycle is present, indicating the existence of bifurcations 
induced by the model parameters. In Figure ga) the 
dissipation constant is increased to £ — 0.02 for which 
a larger limit cycle is observable. Finally, £ w 0.028 is 
the onset of the regime in which all trajectories flow into 
a unique stable cycle and there are no critical points in 
the phase space. This behavior is similar to that of a 
Josephson junction between two superconductors, with 
bias current above the critical current, I = > l,(see 
e.g. Strogatz [E]). There, the current phase 4> and an- 
gular velocity <p = y satisfy the dimensionless equations 



4> = v, 



y = I - sin(0) - ay. 



(14) 



(15) 



where a is the damping parameter. Note that for I > 1, 
there are no critical points. 



A. Codimension-1 Bifurcations 

In order to investigate the bifurcation dynamics, we 
first consider codimension one bifurcations which can be 
induced by varying a single parameter of the model. We 
employ a numerical continuation software called Mat- 
cont |17] , which computes how a critical point advances in 
the phase space by varying model parameters, and char- 
acterizes the point along the continuation according to 
its Jacobian. We begin by the continuation with respect 
to the dissipation constant £■ The fixed parameters are 
kd = 6, A = 0.01, and AE = -0.0025. Figure ga) 
shows the computed continuation curve, on which an 
Andronov-Hopf (H) bifurcation, limit point (LP), and 
neutral saddle point (also labeled as H) are observed. 
These points are characterized with the normal form co- 
efficient in Table [TTJ. The Andronov-Hopf bifucation is 
supercritical, since the first Lyapunov coefficient is neg- 
ative [IHl HH1 IIH]- The vertical lines indicate the span 
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FIG. 3: (Colour online) Phase-space plot of coupled soliton-SP system with population imbalance dissipation. Model parameters arc kd — 6, A — 
0.01, C — 0.01. (a) AE — 0, (b) AE — —0.0025. White curves indicate sample trajectories. The magnitude of the gradient plot is color-coded 
logarithmically. In (b), only trajectories exterior to the limit cycles are shown for clarity. 
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FIG. 4: (Colour onlinc)Phasc-space plot of coupled soliton-SP system with population imbalance dissipation. Model parameters arc kd — 6, A — 
0.01, Ai£ — —0.0025 (a) £ — 0.02. Exterior and interior trajectories shown for left and right limit cycles, respectively, (b) £ — 0.028. White curves 
indicate trajectories. Gradient plot is color-coded logarithmically 



of limit cycles in Z axis, emanating from the Hopf-point. 
Figure[5](b) shows the limit cycles in the Z—<j) — Q plane. 
Interestingly, the limit cycles get bigger in the <f>— Z plane 
for stronger dissipation, before they are terminated by a 
limit point cycle. For reference, the limit cycles in Fig- 
ures [3] and [i] are in the range of Fig. [5jb) . 

TABLE II: Bifurcation points shown in Fig. [5|a) 



Label AE 




Z* 




H 0.0437 


-1.05 


-0.434 


-0.0126 


LP 0.0272 


-1.42 


-0.767 


-0.109 


NS (H) 0.0184 


-1.48 


-0.495 





LP: Limit Point (fold),H:Andronov-Hopf, NS:Neutral 
Saddle, a is the normal form coefficient for LP and first 
Lyapunov coefficient for H. 

Next, the continuation is performed in the asymmetry 
parameter AE. Figures j6^a) and j6^b) show these results 
for £ = 0.0095. In Table |TTT} the bifurcation points are 



characterized and the respective normal form coefficients 
are listed. The first Lyapunov coefficient is negative for 
the Andronov-Hopf points which indicates supercritical 
bifucation. For clarity, the limit cycles are plotted only 
on the AE > side; they are spawn between the two 
Andronov-Hopf points on the AE > side, too. 

In order to gain further insight about the bifurcation, 
we investigate the continuation curve in the Z— AE plane 
for different values of A in Fig. [7] and for £ in Fig. [8] Equi- 
librium points are periodic in the Z — <p plane, thus, each 
continuation curve closes onto itself in these parametric 
planes. The continuation curves extend in a finite region 
of AE — Z plane, where the range also depends on other 
model parameters. 

In Figure [7j we depict the effect of nonlinearity (A) on 
the continuation curve. In view of Eq. [8] we note that, 
for a given soliton amplitude, increasing nonlinearity con- 
fines the soliton more tight around its propagation axis 
and decreases the coupling between the soliton and the 
SP. Thus, while the continuation curve extends in the 
whole Z axis for small A (blue curve) it shrinks and gets 
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FIG. 5: (Colour online)(a) Continuation of the stable point as a function of dissipation constant. Model parameters are kd — 6, A — 0.01, 
AE — —0.0025. Andronov-Hopf (H) bifurcation, limit point (LP), and neutral saddle (H) arc observed. Vertical lines indicate the extend of the 
limit cycles along the Z-axis. A limit point cycle terminates the limit cycles at C, ~ 2.01, (b) Limit cycles plotted in the Z — — £ parametric 
phase- space. 




FIG. 6: (Colour online) (a)Continuation of the stable point as a function of AE. Model parameters are kd — 6, A — 0.01, C — 0.0095. Supercritical 
limit cycles are spawned between the two Hopf-bifurcation points, (b) Visualization of the limit cycles indicated in Fig. [ma) in the Z — <f) ~ A.E 
parametric phase-space. 



TABLE III: Bifurcation points shown in Fig. [6^ a) 



Label 


AE 




Z* 




LP 


-0.340 


0.878 


-0.997 


3.917 


H 


-0.148 


0.984 


-0.946 


-0.00626 


H 


-0.00257 


0.990 


-0.434 


-0.00275 


LP 


0.00319 


-1.01 


-0.139 


0.219 


LP 


0.00909 


-0.311 


0.797 


-0.0148 


H 


0.0112 


0.0102 


-0.4336 


-0.00098 


H 


0.167 


0.0165 


-0.9461 


-0.0623 


LP 


0.360 


0.0102 


-0.434 


-3.92 



LP: Limit Point (fold), H:Hopf point, a is the normal form 
coefficient for LP and first Lyapunov coefficient for H. 



confined to the range close to Z = —1 for strong non- 
linearity (black curve). Note that multiple equilibria are 
present for a given AE value when \AE\ is small (e.g. 
< 0.05). In addition, the nonlinearity also induces dif- 
ferent bifurcation points in the continuation curve. For 



A = 0.003 only LP occur. For A = 0.03 LP and H bifurca- 
tions are present but the first Lyapunov exponent is posi- 
tive (i.e. subcritical Hopf bifurcations) . For the interme- 
diate values between these extremes, we observe that the 
continuation curve hosts multiple H-points which appear 
in pairs in negative and positive sides of AE axis. These 
are supercritical, and as it was shown in Fig. [6] stable 
limit cycles are spawn between them. 

When the dissipation strength £ increases, the contin- 
uation curve in the AE — Z plane shrinks as shown in 
Fig. [5] Limit point and Hopf bifurcations decorate the 
continuation curves. The Hopf points close to Z = — 1 
are neutral saddle points. The other Hopf points are su- 
percritical. As it was shown in Fig. [5] stable limit cycles 
bifurcate from these Hopf points. Although not shown 
here, the curve completely dissappears for £ > 0.028 since 
no critical point is present. 

We recall that AE = A — v p is the difference be- 
tween the (soliton-characteristic) nonlinearity and the 
(SP-characteristic) propagation constant. While not ob- 
vious in these terms, it is analogous to the asymmetry 
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FIG. 7'. (Colour online) AE continuation of the stable point for dif- 
ferent values of A, with kd — 6, Q — 0.005. Hopf (H) and limit point 
bifurcations arc labeled on the continuation curves. 
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FIG. 8: (Colour online) AE continuation of the stable point for dif- 
ferent values of £, with kd — 6, A — 0.01. Hopf (H) and limit point 
bifurcations arc labeled on the continuation curves. Inset shows en- 
largement around the AE — 0. Z — region. Blue curve corresponds 
to the yellow curve in Fig. [7] 



between the wells of the double well potential in BJJ 
systems, and thus, sizes the asymmetry between the soli- 
ton state and the SP state of the Josephson junction. In 
Figures [7] and |5J we observe that for large magnitudes of 
AE (i.e. irrespective of its sign), the equilibrium points 
lie in the SP-dominant ( — 1 < Z << 0) region. Soliton 
dominant equilibrium points occur under weak nonlin- 
earity and weak dissipative conditions only. 



B. Codimension-2 Bifurcations 

Codimension-2 bifurcations require the simultaneous 
variation of two parameters. The blue curve in Fig. [9] 
shows the continuation an equilibrium point in the ^— AE 
parameter space. We observe three Bogdanov-Takens 
(BT) bifurcation points and a cusp point. BT-bifurcation 
implies the existence of three codimension-one bifur- 
cations nearby, namely, a saddle-node bifurcation, an 
Andronov-Hopf bifurcation and a homoclinic bifurcation. 
From the BT points, the continuation curves of the Hopf 
point and limit point can be obtained respectively. The 



orange curve is the Hopf curve passing through the BT 
point close to the cusp point. The green curve is the 
Hopf curve passing through the lower BT point. Each 
Hopf curve has two generalized Hopf-points when they 
intersect the C — axis. The black curve is the continu- 
ation of the limit point from the upper BT point on the 
green Hopf curve. 

This codimension-2 diagram indicates the regions of £ 
and AE in which stable limit cycles can be generated in 
the system. In Figure [6] it was shown that for fixed £ 
stable limit cycles are spawn between two Hopf points. 
Here, this corresponds to the region between the orange 
and the green curves. The Cusp point on the blue curve 
marks the value of ( = 0.028 beyond which there are 
no critical points for A = 0.1 and AE = —0.0025 as it 
was shown in Fig. [5Jb) and FigQb). In Table [Tv] the 
coordinates of the bifurcation points on the blue curve 
are given in the £ — AE plane and <j)~ Z plane and their 
normal form coefficients are listed. BT points have two 
normal form coefficients with respect to two varying pa- 
rameters. Cusp point is characterized by a single normal 
form coefficient. 




FIG. 9: (Colour online) Codimcnsion-2 bifurcation diagram in the 
£ — AE plane. Model parameters are kd — 6, A — 0.01. Blue curve 
shows the continuation of an equilibrium point. Orange and green 
curves arc the Hopf continuation curves from corresponding BT points. 
Black curve is the limit point continuation curve. BT: Bogdanov- 
Takens point, C: Cusp point. GH: Generalized Hopf point. Inset shows 
the zoomed view around the AE — axis. 



V. CONCLUSION 

In this work, we investigated the dynamical features 
of a dissipative nonlinear Josephson junction which is 
formed by a paraxial optical soliton surface-plasmon cou- 
pled system. The most notable part of the heuristic 
model implemented here is the nonlinear nature of the 
coupling mechanism. For a dissipation proportional to 
the angular velocity, the system decays into stable fixed 
points with constant relative phase and constant average 
population imbalance. A particular feature observed here 
is the phase-slip phenomenon, where the odd-7r modes of 
the dissipationless system decay into even-7r modes. For 
dissipation proportional to the population imbalance, su- 
percritical Andronov-Hopf bifurcation can occur which is 
notable in that stable oscillatory modes between the soli- 



TABLE IV: Codimension-2 bifurcations in the £ — AE parameter space (blue curve in Fig M 



Label 


C AE 


4>* 


Z* M) (1) ,c' 2 > 


BT 


3.19 x 1(T 4 


-1.00 


-0.139 (1.51 x 10~ 4 ,-1.46 x 10~ 7 ) 


BT 


0.026 1.16 x 10~ 3 


-1.35 


-0.434 (1.38 x 10~ 4 ,-7.18 x 10" 2 ) 


CP 


0.028 3.61 x 10~ 3 


-1.46 


-0.65 -0.10 


BT 


0.014 -0.096 


-1.27 


-0.946 (-0.12,0.29) 



BT: Bogdanov-Takens Point, CP:Cusp point, 
normal form coefficients of BT. ' 2 ' normal form coefficient of CP. 



ton and surface-plasmon fields can be maintained against 
strong dissipation. The heavily damped limit resembles 
the behavior of Superconducting Josephson junctions in 
which the system decays rapidly into a unique oscillat- 
ing state. In view of this analysis, we propose that the 
nonlinear Josephson junction between the optical soliton 
and surface plasmons can exhibit interesting dynamical 
features which may be exploited further based on more 
rigorous theoretical models. 
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